Landau Fermi Liquid Picture of Spin Density Functional Theory: 
Strutinsky Approach to Quantum Dots 
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We analyze the ground state energy and spin of quantum dots obtained from spin density functional theory 
(SDFT) calculations. First, we introduce a Strutinsky-type approximation, in which quantum interference is 
treated as a correction to a smooth Thomas-Fermi description. For large irregular dots, we find that the second- 
order Strutinsky expressions have an accuracy of about 5 percent compared to the full SDFT and capture all the 
qualitative features. Second, we perform a random matrix theory/random plane wave analysis of the Strutinsky 
SDFT expressions. The results are statistically similar to the SDFT quantum dot statistics. Finally, we note that 
the second-order Strutinsky approximation provides, in essence, a Landau Fermi liquid picture of spin density 
functional theory. For instance, the leading term in the spin channel is simply the familiar exchange constant. 
A direct comparison between SDFT and the perturbation theory derived "universal Hamiltonian" is thus made 
possible. 

PACS numbers: 73.21.La, 73.23.Hk, 05.45.Mt, 71.10.Ay 



I. INTRODUCTION 

Semiconductor quantum dots are now routinely obtained 
using electrostatic gates or etching processes to pattern a two 
dimensional electron gas formed in some heterostructure (typ- 
ically GaAs/AlGaAs).' Many of their sometimes surprising 
properties are now reasonably well understood from a quali- 
tative or statistical point of view,' - and it is now realistic to 
think about using these quantum dots for some specific pur- 
pose, such as spin filtering,^ current or spin pumping,'* or in 
the context of quantum information.^ 

In this context, it becomes important to go beyond a qual- 
itative or statistical description, and to develop tools able to 
predict quantitatively the properties of a specific quantum dot 
for given parameters. For isolated or weakly connected dots 
- the Coulomb Blockade regime - the Coulomb interaction 
between electrons plays an important role and has to be taken 
into account properly. 

A method of choice in this regard is the density functional 
approach,*"'^ which has been widely used in theoretical model- 
ing of quantum dots.*^ More specifically, since it is necessary 
to describe correctly the spin degree of freedom, we have in 
mind a spin density functional, where each density of spin 
n°'{r), with a = a,i3 corresponding to majority and minority 
spins, are treated as an independent variable. How this can be 
achieved in practice for a dot containing up to four hundred 
electrons, and for an parameter as high as four, has been 
demonstrated in a series of papers.^ '"'" 

One striking feature of these calculations, however, is that 
the qualitative picture which emerges is somewhat unex- 
pected in view of previous results. Indeed, within a statis- 
tical approach and assuming the classical dynamics within 
the nanostructure is sufficiently chaotic, one can model the 
wave functions in the quantum dot using Random Matrix The- 
ory (RMT). If furthermore the Coulomb interaction is treated 
within the Random Phase Approximation (RPA), it is possible 
to derive various statistical quantitiesjiSi^iiiii^ such as the dis- 



tribution of spacing between Coulomb Blockade conductance 
peak, or the probability of occurrence of non-standard spins 
[that is, not zero (not one-half) for even (odd) particle num- 
ber]. It turns out, for instance, that the spin density functional 
calculations give a larger number of "high" spins than was 
predicted within RMT plus RPA modeling.'' Such discrepan- 
cies could originate from a variety of causes, ranging from the 
the statistical behavior of the Kohn-Sham wave functions to 
an intrinsic failure of one or the other approach. The goal of 
this paper is to clarify this issue. 

For this purpose, we need a way to understand, or to orga- 
nize, the numbers obtained from the full-fledged spin density 
functional theory (SDFT). This can be done using the Strutin- 
sky approximation to DFT discussed in Ref. 16, up to straight- 
forward modifications to include the spin variable. 

The Strutinsky approach is an approximation scheme, de- 
veloped in the late sixties in the context of nuclear physics, 
in which the interference (or shell) effects are introduced 
perturbatively.i^ It has been since then used in many sub- 
field of physicSfi^ including the calculation of the binding 
energy of metal clusters.^" Standard treatments are usually 
limited to first-order corrections; in this case, it is similar to 
the HaiTis functional^' familiar to the DFT literature. Here, 
in contrast, we use the second-order approximation^^ devel- 
oped in Ref. 16. We shall see that the second-order Strutinsky 
scheme turns out to be extremely accurate in some circum- 
stances. Furthermore, even when it is less precise, which hap- 
pens in conjunction with the occurrence of "spin contamina- 
tion" in the SDFT calculations, it still provides a qualitatively 
correct statistical description. 

The second-order Strutinsky corrections can be cast in a 
very natural form: '^ in fact, they amount to taking into account 
the residual (screened) interactions between quasiparticles in 
a Landau Fermi liquid picture. They are thus amenable to 
treatment by the same RMT approach as was used previously 
for RPA, providing an analogue of the "universal Hamilto- 
nian" in the case of SDFT. Since within the Strutinsky approx- 
imation the quantum dot properties are relatively transparent. 
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we will then be in position to discuss the difference between 
SDFT results and those obtained from RMT plus RPA. 

The paper is organized as follows. In Section|ll| we briefly 
review the Strutinsky approximation as it applies to spin den- 
sity functional theory. In Section|lIlJ we consider more specif- 
ically the local density approximation from this perspective, 
and in particular the screened potential that it implies. Sec- 
tion II VI covers in detail the specific case of a model quan- 
tum dot with quartic external confining potential. This model 
is used to investigate the accuracy of the Strutinsky approxi- 
mation for various electron densities. In Section |V] we dis- 
cuss how the Fermi liquid picture emerging from the Strutin- 
sky approximation scheme can be used, in conjunction with 
a random plane wave model of the wavefunctions, to ana- 
lyze the peak spacing and spin the distributions resulting from 
the SDFT calculation. This framework also makes it possible 
to discuss the mechanism of spin contamination. Finally, in 
the last section we come back to the original question moti- 
vating this work and make use of what we understand from 
the Strutinsky approximation to discuss the discrepancies be- 
tween SDFT calculations and RPA plus RMT predictions. 



n. STRUTINSKY APPROXIMATION FOR SPIN DENSITY 
FUNCTIONAL THEORY 

In the spin density functional description of a quan- 
tum dot, one considers a functional of both spin densities 

[?i"(r),n'3(r)] 

^KsK,n'3] =TKsK,n'5]+£:t„tK,n'5] , (1) 

where a = a,P correspond to majority and minority spins, 
respectively. In this expression, the second term is an explicit 
functional of the densities. 



£tot[n",n'^] = Jdr n{r) C/c.t(r) 

+ / drdr' n(r) «int(r - r') n(r') + n''] , (2) 



where n(r) = n"(r) + n^{r), J7cxt(r) is the exterior con- 
fining potential, and the precise form of the exchange corre- 
lation term Sy^dn" , n^] is to be discussed in more detail in 
Section|lIl| tiint(r,r') is the electron-electron interaction ker- 
nel. The presence of metallic gates can be taken into account 
by including an image term in addition to the bare Coulomb 
interaction. 



(3) 



where Zd is the distance between the top confining gate and 
two-dimensional electron gas. The image term in the inter- 
action kernel greatly reduces classical Coulomb repulsion be- 
tween electrons at a distance larger than Zd so that the elec- 
tron density far from the boundary is quite uniform. The bare 
Coulomb interaction can be recovered by letting Zd go to in- 
finity. 



The kinetic energy term Tusin" , n^], on the other hand, 
is expressed in terms of the auxiliary set of orthonormal 

functions vl"'''^ (* — 1, • ■ ■ i -^a,/?) such that n°'(r) = 



Et'il^fWPas 

•' a=a,f3 i=l 

^From the density functional Eq. Q , the ground state en- 
ergy of the quantum dot containing {Na, Np) particles of spin 
(a, (3) is obtained as 



(5) 



where the Kohn Sham densities [n'^^{v)^n^^{v)] minimize 
J-yis under the constraint given by the number of particles of 
each spin. This in practice is equivalent to solving the Kohn- 
Sham equations 



-— + C/ks (r) ) = e^CW , ^ = 1, -.N, 



with the spin dependent self-consistent potential 



(5n'^(r) 



tot r a 

IIkC, I It, 



'■KSi "KSJ 



(6) 



(7) 



In this section, we shall give a brief description of the 
second-order Strutinsky approximation as it applies to spin 
density functional theory. Up to the introduction of the spin 
indices, the derivation of this approximation follows exactly 
the same lines as the spinless case discussed in details in Ref. 
Ilfil We shall therefore not reproduce it here, but rather try to 
emphasize what exactly are the assumptions made in deriving 
the approximation, and then merely write down the expression 
we shall use in the following sections. 



A. Generalized Thomas-Fermi approximation 

The basic idea of the Strutinsky energy correction method 
is to start from smooth approximations, Etf and ?T.TF(r), to 
the DFT energy Eks and electronic density riKs(r). Then, 
fluctuating terms are added perturbatively as an expansion in 
the small parameter 



i{r) = nKs(r) - '^TF(r) 



(8) 



In the original work of Strutinskysi^ii^ various ways of con- 
structing the smooth approximation have been considered. 
The most natural choice for n^p{r) turns out to be the so- 
lution of the Thomas-Fermi equation 



(9) 



coupled with £^TF = •^tf['t-tf('")' "'tf(^)]- The "general- 



ized" Thomas-Fermi functional is defined as 



(10) 
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with Stot [n", n^] given by Eq. It differs from the original 
spin density functional only in that the quantum mechanical 
kinetic energy Tks, Eq. ©, is replaced by an explicit func- 
tional of the density. For two dimensional systems this takes 

the form Ttf ,n'^]= T^p ' "'^] + "^tf ' "'^1 with 



1 



2N{0) 
87rA^(0) 



dr[n"(r)2 + n'3(r)2] (11) 



dr 



12) 



where N{0) ~ m/Trh? the density of states for 2-dimensional 
systems and 77 a dimensionless constant taken here to be 0.25. 
Such a choice for the kinetic energy functional correctly takes 
into account the Pauli exclusion principle, and thus that an in- 
crease in kinetic energy is required to put many particle at the 
same space location, but fails to include fluctuations associ- 
ated with quantum interference. 

To illustrate this, let us for a short while the assume 77 = 
0, i.e. only the first term T^p of the Thomas-Fermi kinetic 
energy is taken into account. Then, one can show that the 
Thomas-Fermi density fulfills the self-consistent equation 

n^F(r) = n^[[/-p](r) (13) 

where the Thomas-Fermi potential is defined, as in Eq. Q, by 



C/^F(r) 



(5n'^(r) 



and 



raw 



dp 



(27ra)' 



e 



P 

2m 



t/?F(r) 



(14) 



(15) 



is the Weyl part of the density of states of a system of inde- 
pendent particles evolving under the potential L'^p(r) (0 is 
the Heaviside step function, d = 2 is the dimensionality, and 
the chemical potentials /i^p are chosen such as to fulfill the 
constraints on the total number of particles with spin a and 
13). From its definition, n°'[C/fp](r) (and as a consequence 
n^p(r)) is a smooth function, in the sense that it can change 
appreciably only on the scale on which ?7cxt(r) varies, but 
cannot account for the quantum fluctuations of the density on 
the scale of the Fermi wavelength. 

The Weyl approximation is, however, only the leading term 
in a semiclassical expansion of the smooth part of the density 
of states, and higher-order corrections in h can be added in a 
systematic way. The standard way to choose T^p is such 
that Eq. M3\ holds, but with an approximation to the smooth 
density of states fi°'[f/,Jp](r) which includes both the Weyl 
term Eq. (I15> and the first h corrections. For two dimensional 
systems, however, the corrective term computed from this pre- 
scription turns out to be zero, when the presence of 7^p is 
actually useful in smoothing the Thomas-Fermi density near 
the boundaries of the classically allowed region. We have 
therefore used the phenomenological Weisacker-like terrain 
Eq. (I12t which plays a similar role. 



B. Strutinsky correction terms 

The practical implementation for the Strutinsky scheme can 
be summarized as follows. The first step consists in solving 
the generalized Thomas-Fermi equation Eq. (|9}, which de- 
fines a zeroth-order approximation for the ground state energy, 
E'YY, as well as an approximation to the density of particles 
nf.p(r). From this latter quantity, one derives for each spin 
a = a,(3 the effective potential J7fp(r) through Eq. MAI . 

The second step consists in solving the Schrodinger equa- 
tions (again for each spin) 



2m 



V2 + C/^p(r) 



^^f(r), i = l,...,N^. 

(16) 

Therefore, while the Kohn-Sham equations are both quan- 
tum mechanical in nature and self-consistent, here all the self- 
consistency is left at the "classical-like" level of the Thomas- 
Fermi equation, and the quantum mechanical wave interfer- 
ence aspect is taken into account without self-consistency. 
One obtains in this way a new density 



(17) 



It can be seen'^ that n!^p(r) is by construction a smooth ap- 
proximation to n"^ (r ) [this is basically the content of Eq. ( I13H . 
Therefore 



<,(r)EEn-(r)-n^p(r) 



(18) 



can be considered as the oscillating part of n"', and will de- 
scribe the short scale variations of the density associated with 
quantum interference. 

Once the eigenfunctions and eigenvalues of Eq. il6\ are 
known, corrections to the Thomas-Fermi energy can be added 
perturbatively 



Eks - Etf + '^E^^^ + AE^^^ 



(19) 



The first-order correction turns out to be simply the oscillating 
part of the one particle eneigyi^ii^ 



with 



Ni3 



?/3 



(20) 



(21) 



i=l 



the one particle energy, and 

fip ^TtfKf.^f] + E JdrU!}p{r)n^p{r) (22) 

its smooth part. 

The second-order correction AE''^' can be expressed in two 
different waysi^ 

^E'^'^^l E jdvdv'fil^,{v)V^£{vy)6n^\v') 



(23) 
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(7",cr'— a,/3 

(24) 

with 



T rcr.cr / i\ ^ ^tot r a /3 1 /'^C\ 

Kare(r,r ) = fa.(,)fa.-(,,) ["TF,^TF] (25) 



(T,CT 

screened 



(r,r') 



E 

■"=a,/3 • 



TF 



Vba 



(r,r") 



iTF 



T4a 



(r",r') 



(26) 



(The matrix inversion here should be taken with respect to 
both the spatial position and the spin indices.) Note that in 
Eq. M6\ we can use 



(0) 



TF 



(5n'^(r)(5n<^'(r') A^(0) 



Sa.a' S{r' - r) 



(27) 



and, neglecting terms involving derivatives of the Thomas- 
Fermi density, 



(1) 



TF 



V 



■Ar'S{r' 



Sn<^{r)6n'^'{r') 47r7V(0) n^p(r) ' ^^^^ 

The first of these second-order expressions, Eq. ( I23t . in- 
volves the 5n'^{r) which are in principle unknown. It is there- 
fore not useful if one is trying to avoid performing the full 
Kohn-Sham calculation. Here, however, we shall also do this 
full self-consistent calculation in order to evaluate the accu- 
racy of the Strutinsky scheme; we therefore know Sn"'{r). In 
this case Eq. ( I23> turns out to be somewhat simpler to com- 
pute than Eq. i24\ . This is not, however, because of the ma- 
trix inversion in Eq. i26\ - for the case we shall consider, the 
screening length is much smaller than the size of the system, 
allowing this inversion to be done analytically. Rather, the re- 
sulting screened potential will turn out to be diagonal neither 
in the position nor in the momentum representation (see next 
section), in contrast to Vb'aro which is the sum of two parts, 
one diagonal in the position and the other in the momentum 
representation. Since Eq. ( I23> and ( I24> are essentially at the 
same level of precision, we shall in the following mainly use 
the first one Eq. (I23> : we will apply Eq. (I24> only with an ap- 
proximate 14crconcd which is diagonal in the momentum rep- 
resentation, thus inducing additional errors that are not intrin- 
sic to the Strutinsky scheme. 

C. Applicability of second-order Strutinsky expressions 

Although we shall not reproduce here the derivation of 
Eq. (I23> (see again Ref.Hfilfor details), the condition of appli- 
cability of this equation can be understood easily by compar- 
ing the Kohn-Sham equations Eq. (|6j with the ones defined by 



the Thomas-Fermi potential U^p{r), Eq. il6\ . We see there 
that the Kohn-Sham wave functions ipf and their approxima- 
tions are defined through Schrodinger equations that differ 
only by a difference in the potential term 



(29) 



E 

''=a,l3 



dr'V:^:Ar,r')6n'^ (r') 



with Vbarc(r, r') defined by Eq. i25\ . 

The main approximation in the derivation of Eq. i23\ is that 
the SU (r) can be taken into account by a second-order pertur- 
bative calculation. In general, this implies that the non diag- 
onal matrix elements {5U)fj = (^f |5t7'^|0J) {i ^ j) should 
typically be smaller than the mean level spacing A between 
one particle energies. More specifically, since the only rele- 
vant errors are ones modifying the Slater determinant formed 
by the occupied orbitals, the actual condition is that the matrix 
element between the first unoccupied orbital and the last oc- 
cupied orbital is smaller than the spacing between these level, 
that is 

{rN^+i\SU"\rNj « (ejv^+i-ejvj = a,(3). (30) 

A good accuracy, on the scale of the one particle mean level 
spacing, of the Strutinsky approximation in the form Eq. (I23> 
is equivalent to the above condition. 

Now, Eq. ( I30> again involves 5n°'{r), which is unknown. It 
is therefore not possible to prove rigorously that it should be 
fulfilled. It is possible, however, to show the consistency of 
such an assumption. If we assume that 

(i) |(0f |(5J7'"|0J)| < A for all pairs of orbitals and 
a = a,(3, so that 6U"' can be treated perturbatively, and 

(ii) the oscillating parts of the Kohn-Sham densities 
('^Ks)osc(i") do not differ significantly from ^^^^(r), 
implying that Sn"'{r) ~ 6n"'{r) + n'^^^{r), 

then it can be shown that 



(a) 



a' —a,j3 



dv' 



ed(r,r')nosc(r') (31) 



which immediately yields Eq. J24> from i23\ . and 

(b) (<5t/2.)/A2 and ([(nKs)^.e(r) - <.c(r)]^)/([5n-(r)]^) 
are both negligible for a large dot for which the chaotic 
motion in the potential U^p allows modeling of the 
wave functions in terms of random plane waves (being 
of order \ng/g with g the dimensionless conductance 
of the dot). 



III. THE SCREENED INTERACTION IN THE LOCAL 
DENSITY APPROXIMATION 

The second order Strutinky corrections, Eqs. ( I23> or J24t . 
involve the bare and screened interactions Eqs. (I25> and i26\ . 
In this section, we shall consider the particular form these in- 
teractions take for the exchange correlation term we use to 
perform the actual SDFT calculations, namely the local spin 
density approximation 

£^c[n",n^] ~ jdr n{r) e^,{n{r)X{r)) , (32) 

where ({r) = (n" — n^)/n is the polarization of the elec- 
tron gas, and e^c is the exchange plus correlation energy 
per electron for the uniform electron gas with polarization 
We furthermore use Tanatar and Ceperley's form of Cxc at 

C = and IjS, and the interpolation e^c{n, C) = ^xc{n, 0) + 
/(C)[^xc('T-, 1) — exc("-,0)] for arbitrary polarization, with 
/(() = [(1 + C)3/2 + (1 _ ^)3/2 _ 2] /(23/2 _ 2). Recently, 

Attaccalite et ali2i parameterized a new LSDA exchange- 
correlation functional form based on more accurate quantum 
Monte Carlo calculations (see e.g. Ref. that include spin 
polarization explicitly. We have, however, checked that for the 
quantities in which we are interested here, this new functional 
introduces only minor modifications; we shall therefore use 
Tanatar-Ceperley's form for our discussion. 

From the expression of the functional E^dn", n^], the bare 
and screened interaction potentials Eq. ( I25> and(l26>. needed 
for the evaluation of the second-order Strutinsky corrections, 
are easily computed. The bare interaction is the sum of two 
terms Vbarc — V^oui + 14c- The Coulomb interaction 



1.5 



Fcoui(r,r') =Vi„t(|r-r'|) 



1 1 
1 1 



(33) 



is independent of both the density and spin. The matrix struc- 
ture here refers to the spin indices a and /3. The exchange 
correlation term is local and can be expressed as 



14c(r,r') = -27raoe25(r-r') 

va{n{r),ar)) Vb{{n{r),ar)) 
Vb{n{r),C{r)) Va{n{r), -({r)) 



(34) 



where Va and Vb are obtained from the partial derivative of 
exc(f^, C) [defined in Eq. ( I32H . In all numerical calcula- 
tions we shall keep entirely the dependence of Va and Vb on 
the polarization C,. However, this latter will usually be rela- 
tively small, and Va and Vb contains already second deriva- 
tives of the functional £xc with respect to the polarization. 



— v^(r,C=0) 

- - v,(r ,C=0) 




FIG. 1: (Color online) The functions Va{n,0) and Vb{n,0), which 
define the functional second derivative of E^c [n" , n*^] [see Eq. <34H . 
as a function of the parameter rs = 1/ 110^71. 



We shall therefore proceed assuming Va{n, C) — Vain, 0) and 
Vb{n,C) — Vb{n,0). The dependence of these functions on 
the parameter — (vraQn)^^/^, with ao — It? /me^ the 2D 
Bohr radius, is shown in Fig^ 

Turning now to the screened interaction K!creenod(r, r')' 
is useful to switch to the variables (R, 1) = ((r' + r)/2, (r — 
r')). The Fourier transform of Vbarc (R-, 1) with respect to 1 
reads 



Vbare(R,q) = 27re^aot;c(q) 



1 1 
1 1 



(35) 



27raoe 



Vb{K) w,(R) 



where 



ao|q| 



(36) 



(again, the pure Coulomb case is recovered by letting the dis- 
tance to the top gate Zd go to infinity). This can by further 
simplified by diagonalizing the matrix in the spin indices: the 
eigenvectors are the "charge channel" c = (a + \/2 and 
the "spin channel" s — [a — /3)/V2. The eigenvalues are 



Ar-(R,q) 



N{Q) 

in the charge channel and 
A^-(R,q) = 



{2v,{q) - [Va(R) + Vbi'R)]} (37) 



iV(0) 



K(R)-Va(R) 



(38) 



in the spin channel [note 27raoe^ ~ 2/A^(0)]. 

Because the screening length is short on the scale for which 
the smooth part of the electronic density varies appreciably, 
the operator S'^Tt^p/Sh? + Vbarc can be inverted by treating 
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the variable R as a parameter, appearing thus diagonal in the 
q representation. One obtains in this way 



scrccnG( 



i^^'^^- l^(A,-A,)/2 (A, + A,)/2 ) 



(39) 



in terms of the eigenvalues 

2 2w,((z)-K(R) + Wf,(R)] 



Ae(R,q) 



7V(0) 1 + g{q)[2v,{q) - («,(R) + Vb(Ii))] 

(40) 

for the charge channel and 



A.(R, q) 



K(R)-i;a(R) 



7V(0)1+5(q)K(R)-«,(R)] 



(41) 



for the spin channel. The R dependence of Ac and As arises 
from Va{n{IV)) and Wb(ri(R)), and therefore from the local 
value of the density (that is, of the parameter rg) at the location 
the interaction is taking place. Furthermore, the function 



g{q) =1 + 77 



87rnTF(R-) 



(42) 



would just be 1 in the absence of the h correction 7^p to 
the Thomas-Fermi kinetic energy term; it prevents effective 
screening to take place for large momenta. 



IV. THE GATED QUARTIC OSCILLATOR MODEL 

To evaluate the accuracy of the Strutinsky approximation 
scheme, we consider a two dimensional model system for 
which the electrons are confined by a quartic potential 



y'^x)r 



(43) 

The role of the various parameters of J7cxt(a;, y) is the follow- 
ing: a controls the total electronic density (i.e. the parameter 
Ts), and therefore the relative strength of the Coulomb inter- 
action; A allows one to place the system in a regime where the 
corresponding classical motion is essentially chaotic; finally, 
6 and 7 have been introduced to lower the symmetry of the 
system. In the following sections, we use [al] to denote the 
parameter value a = 10^^ x En/a^ (with En = e^/2ao), 
and [a4] for a = IQ-^ x Er/o^, which at iV = 120-200 cor- 
respond to Ts — 0.3 and 1.3, respectively. We use A = 0.53 
and 7 — 0.2 in our calculations unless specified otherwise. In 
addition to this 2-dimensional potential, we assume the exis- 
tence of a metallic gate some distance Zd away from the 2D 
electron gas whose purpose is to cut off the long-range part of 
the Coulomb interaction. This gate is placed sufficiently far 
from the electrons to prevent the formation of a density deficit 
in the center of the potential well without modifying qualita- 
tively the quantum fluctuation. In practice, we take Zd about 
0.75 ao for [al] , and 2.5 for [a4],^ 



A. Electronic densities 

For any set of parameters defining the potential Eq. J43t . 
and for any number of up and down electrons [N°',N^), 
we can compute the Kohn-Sham energies E]^^[N" ^ N^] and 
densities n^^fr) following the approach described in de- 
tail in Ref. 10. For the Strutinsky approximation, the only 
part which presents some degree of difficulty is actually 
the Thomas-Fermi calculation, for which we have devel- 
oped a new conjugate-gradient method which turns out to 
be extremely efficient3§^ Once ni^'^'' (r) are known, the ef- 
fective potentials U^^\v) and the corresponding densities 
h^°''^\r) follow immediately. 

Figure |2] shows the densities nTF(r), nKs(r), (5n(r), and 
riosc (r) for the ground state with N = 200 electrons (Na — 
Np = 100) of the gated quartic oscillator system with pa- 
rameter [a4], corresponding to an interaction parameter of 
Ts — 1.3. Noting in particular the difference of scale between 
the upper and lower panels, one can observe that the Thomas- 
Fermi density already is a very good approximation to the ex- 
act one, and that 6n{r) is a small oscillating correction, of the 
same magnitude as the oscillating part of h{r). Apparent also 
on this figure is the fact that the largest errors are located at the 
boundary of the dot where corrections to the Weyl density of 
higher order in Ti are the largest. To make this more visible, we 
plot in Fig.l^the densities dn{r) and nosc(r) along a cut on a 
diagonal of the confining potential for two sets of parameters. 
This makes it clear that rjosc (r) is an oscillating function only 
in the interior of the dot, but has a proportionally large secular 
component at the boundary. As a consequence, choosing cor- 
rectly the term T^p of the Thomas-Fermi kinetic energy term 
is actually important to obtain good accuracy. We have there- 
fore determined the parameter rj — 0.25 of this functional by 
imposing that 8°^^ [N] oscillates around zero, rather than hav- 
ing a significant mean value. As an illustration, we also plot in 
Fig.|3lthe same quantities but for a calculation where a value 
?; = 1 has been used for the Thomas-Fermi kinetic energy 
correcting term. We see that this increases significantly the 
error in the Thomas-Fermi density at the boundary, reducing 
the effectiveness of the Strutinsky approximation. 



B. Ground state energies 

With the expressions Eqs. i35i and ( I39I I of the bare and 
screened interactions, and using the known eigenvalues and 
eigenfunctions for the Schrodinger equations Eqs. (|6jl and 
( I16t . the Strutinsky approximation for the total energy, includ- 
ing the order one Eq. i2Q\ and order two Eqs. ( I23t or i24l cor- 
rections, can be computed for any {Na, N/s). As for the full 
spin density calculations, the ground state energy for a given 
total number of electron N is then obtain as the minimum over 
{Na,Nfj) with Na + N[3 = N of these energies. 

Let us now consider these ground state energies for a choice 
of parameters [al], such that the coefficient = 0.3 is still 
smaller than one and thus the Coulomb interaction is not 
very large compared to the kinetic energy of the particles. 
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FIG. 2: (Color online) Particle density for the parameters [a4] for a 
system of iV = 200 electrons (Na = Np = 100). From top to 
bottom: nTF(r); nKs(r); 5n(r); and ?7,osc(r). Note that jitp is a 
smooth approximation to nKS and that 5n{v) and nosc(r) are very 
similar. 
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FIG. 3: (Color online) Diagonal cut of particle densities 5n[v) (solid) 
and nosc(r) (dashed), in units of the average Thomas-Fermi den- 
sity riTF inside the dot. Upper panel: parameters [al] (r^ — 0.3); 
Lower panel parameters [a4] (r^ ~ 1.3). The dot contains A'^ = 200 
electrons {Na ~ N/s = 100). The thin dashed lines correspond 
to nose (r) from a less accurate choice of the Thomas-Fermi kinetic 
energy term, namely rj = 1.0 instead of 0.25. 



Fig. |4] displays the corresponding difference between Struti- 
nsky and Kohn-Sham ground-state energies, in mean level 
spacing units, for a number of electrons ranging from 50 to 
200. For the upper panel. Est [N] is obtained from the ex- 
pression Eq. ( I23> using the bare interaction and requiring the 
knowledge of the exact 6n'-°'-^^ . For the lower panel, Est' [N] 
is obtained from Eq. i23\ . which involves the approximate 
screened interaction, but only the knowledge of nose'''' (r). 

The first observation that can be made on that figure is that 
the second form of the Strutinsky approximation appears sub- 
stantially less accurate than the first one. As mentioned ear- 
lier, this is, however, probably due to the fact that, because 
our code was devised to handle only two body interactions 
that were diagonal either in position or in momentum repre- 
sentation, we had to use in that calculation an approximation 
where, for the screened interaction, the local value of the pa- 
rameter Ts (r) = -^/l/TrnTF (r) has been replaced by its mean 
value . 

Indeed, a second feature immediately visible on Fig.l^is 
the presence of a net trend in the energy differences between 
Kohn-Sham and Strutinsky calculations. This secular term is 
related to the non-oscillating component of nosc(r) and 6n{r) 
visible on the lower panels of Figs. |2]and|3]at the boundary 
of the quantum dots. This can be checked by using a less ac- 
curate Thomas-Fermi approximation (e.g. with rj — 1 for the 




Electron Number 

FIG. 4: (Color online) Difference between Strutinsky and Kohn- 
Sham ground state energies, in units of the mean level spacing, as 
a function of the number of electrons, for the high density dot [al] 
(rs — 0.3). Upper panel: Est[N] obtained from Eq. <23> : Lower 
panel: Est* [N] obtained from Eq. <24t . In both cases, the differ- 
ence shows small fluctuations (a few percent) about a linear secular 
trend (dashed line is fit). 

correcting term of the kinetic energy T^p ), and noticing that 
this secular term in the deviation increases noticeably while 
the fluctuations remain less affected. Since the total density 
of electrons at the boundary of the dot is significantly lower 
than its average value, it is relatively natural that in our ap- 
proach, this secular deviation is made significantly worse in 
the second form of our approximation. 

The secular deviation is not completely negligible in terms 
of the mean level spacing. The relevant scale for the smooth 
part of the ground energies is, however, the charging energy, 
compared to which these secular corrections are extremely 
small. To focus on the fluctuating part, we therefore remove 
the secular term (i.e. the straight lines in Fig.|4}, obtaining in 
this way Fig.|5] 

Let us consider first the upper panel, where the first form of 
the Strutinsky approximation has been used. We observe that 
the fluctuating part of the error is usually extremely small, 
typically of order a percent of a mean level spacing, and that 
this is mainly an oscillation between odd and even number of 
particles in the system. Nevertheless, in a few circumstances 
significantly larger deviations are observed, with a magnitude 
typically five percent of a mean level spacing and a sign which 
is always negative. 

To understand the origin of these larger deviations, it is use- 
ful to correlate them with the occurrence of spin contamina- 
tion, that is to situations where the SDFT calculations break 
the spin rotation symmetry. Since the actual spin is a some- 
what ill-defined quantity in a spin density calculation, we need 
however first to specify what we understand by this. Indeed, 
in spin density functional theory the difference A^^ — Np can 
be interpreted as twice the component Sz of the quantum dot 
total spin. Another quantity that can be easily computed is the 
mean value S'(S' + 1) of the operator 5^ for the Slater determi- 
nant formed by the Kohn-Sham orbitals t/jf, i = 1, • • • , N^-, 
which can be expressed as S{S + 1) = Sz{Sz + 1) + 5S with 
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Electron Number 

FIG. 5: (Color online) Difference between Strutinsky and Kohn- 
Sham ground state energies as in Fig. |4] but now with the secular 
trend removed. Solid: all ground states; Dots: ground states with- 
out significant spin contamination. Note the excellent agreement in 
the case of the first Strutinsky form Est (upper panel) - of order 1 
percent - when spin contamination is not present. In the case of the 
second form Est* , the agreement is still very good (lower panel). 

the "spin contamination" 6S given bjsS 

ss^Np- i(^ri^f)r (44) 

From this expression, one sees that if all occupied /3-orbitals 
are identical to the corresponding a-orbitals, SS = 0, and it is 
presumably reasonable to interpret {Na ~ ^/3)/2 as the sys- 
tem's total spin. However, as soon as different-spin orbitals 
start to differ, 6S can take any positive value smaller than Np, 
signaling that, at least, there is some ambiguity in the assess- 
ment of the total spin of the system. In Fig.|6]we have plotted 
the ground state spin contamination (55 [iV] as a function of 
the particle number For [al], the spin contamination is usu- 
ally negligible, except in some few cases where SS's of order 
one half or so are encountered. 

Coming back to the Strutinsky approximation, we can use 
the information from Fig.|5]to exclude in Fig. |5] the ground 
states with significant spin contamination. The remaining 
points correspond to the dotted symbols in this figure. In the 
upper panel, we see that there is almost a one to one corre- 
spondence between larger errors and spin contamination. 

Turning to the lower panel in Fig.|5] we see that the further 
approximations in treating the screening used in evaluating 
Eq. i23i do degrade the accuracy of the ground state energy 
somewhat. Still, ST* gives the fluctuating part of the energy 
to within a few percent. For the spin contamination, no partic- 
ular correlation is seen, presumably again because the overall 
agreement is slightly spoiled by the approximation we made 
for the screened Coulomb interaction. 

In lower density (larger Vg) more realistic dots modeled by 
the parameter set [a4], spin contamination in KS ground states 
is much more pronounced, as shown in the lower panel of 
Fig.|6]- it is, in fact, always significant. In conjunction, Fig.0 
shows that the accuracy of both ST [i.e. using Eq. ( I20l il and 
ST* [i.e. using Eq. ( I23H becomes worse. As at higher density. 
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FIG. 6: Spin contamination 5S[N] as a function of tiie number of 
particles. The spin contamination is infrequent at high density ([al], 
upper panel), but becomes frequent and substantial for rs — 1.3 
([a4], lower panel). 



the main error is a secular trend: in the case of ST* it attains 
a value of several mean level spacings, due presumably to the 
approximations made in treating the screening. After remov- 
ing the secular deviation, however, the fluctuation in the errors 
of ST and ST* ground state energy is still quite small: the rms 
is 0.05A in ST and 0.06A in ST*. Thus for characterizing the 
fluctuating part of the ground state energy, ST* is nearly as 
good as ST, a property we expect to remain valid at larger r^. 



C. Coulomb Blockade peak spacings and spin distribution 

In the previous subsection, we have considered the accu- 
racy of the Strutinsky approximation for individual ground 
state energies. We found that as long as no significant amount 
of spin contamination is present in the SDFT calculation, the 
Strutinsky result provides an excellent approximation when 
Eq. ( I23> is used, and a good, though slightly degraded one, 
with Eq. i24\ . In this latter case, it is probable that the addi- 
tional errors come mainly from the neglect of the local den- 
sity dependence in the screened interaction rather than to the 
Strutinsky approximation itself. We shall come back to the 
issue of spin contamination in the next section, and sketch an 
extension of the theory that would make it suitable to deal with 
the spin contaminated case. 

Before doing so, however, we shall address another ques- 
tion, namely how well, even in the case where a one-to-one 
comparison of ground state energies can imply an error of a 
fraction of mean level spacing, are the statistical properties of 
the quantum dots described within the Strutinsky approxima- 
tion. For instance, we have in mind the distribution of ground 
state spin 5*^ [N] or of ground state energy second difference 
s[N] = E^s[N + 1] + E^s[N - 1] - 2E^s[N + 1]. This 
latter quantity is accessible experimentally by measuring the 
spacing between conductance peaks in the Coulomb Blockade 
regime, and will be referred to below as the "peak spacing". 

In Fig.|8] both peak spacing and spin distributions are plot- 



FIG. 7: (Color online) Difference between Strutinsky and Kohn- 
Sham ground state energies, in units of the mean level spacing, as 
a function of the number of electrons, for the rs ~ 1.3 dot [a4]. 
Upper panel: iSsxfA'']; Lower panel: Est* [N\. The linear secular 
trend (dashed line is a fit) is now significantly larger than for [al] 
(compare with Fig. but the fluctuation about this trend remains 
small (of order 5 percent). 



ted for two interaction strength regime ([al] and [a4]) using ei- 
ther Kohn-Sham results or one or the other form of the Struti- 
nsky approximation. In the small rg case, the agreement is 
naturally excellent, but we see that even for the higher rg case, 
both forms of the Strutinsky approximation give a fairly good 
approximation - certainly they provide a qualitatively correct 
description. 



V. FERMI LIQUID PICTURE 

It is not possible to develop a real Landau Fermi liquid the- 
ory for quantum dots because the mesoscopic fluctuations pre- 
vent any Taylor expansion of the free energy in terms of oc- 
cupation number (Landau theory basically assumes that the 
excitation energies are the smallest energy scales of the prob- 
lem, which is clearly not true for mesoscopic systems because 
of variation on the scale of the mean level spacing.) How- 
ever, discussing what we may call a Landau Fermi liquid "pic- 
ture", in the sense that the low energy physics is described by a 
renormalized weak interaction, is still something meaningful. 

In that sense, what Eqs. (I19>-( I26> express is that SDFT, in 
the limit where the Strutinsky approximation scheme is ac- 
curate, is equivalent to a Landau Fermi liquid picture, where 
quasiparticles with spin a — a^j3 evolve in the effective po- 
tential U^-p and interact through a residual weak interaction 



screened 



(r, r') that can be taken into account as a perturba- 
tion. The only unusual feature is the absence of an exchange- 
like contribution to the total energy. Indeed the main role of 
the exchange correlation functional fxci'T-a, np] is to make the 
interaction between same spin particles different from the one 
between opposite spins. 

Since moreover we have chosen the confining potential 
in such a way that the classical motion within our model 
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FIG. 8: (Color online) Spin and peak spacing distributions for the 
cases [al] (left column) and [a4] (right column). Solid: even A^; 
Dashed: odd N. The statistics are obtained for iV = 120-200 
with (A, 7) = (0.53, 0.2), (0.565, 0.2), (0.6, 0.1), (0.635, 0.15) and 
(0.67, 0.1). From top to bottom: Kohn-Sham , first, and second form 
of the Strutinsky approximation. Agreement between the three meth- 
ods is excellent, of course, for [al]. But even for [a4] where individ- 
ual energies may be in error, the agreement of the distributions is 
very good. 



quantum dot is chaotic, we know that we can use a statisti- 
cal description of the eigenlevels and eigenstates of i?TF = 
P^/2to + J7TF(r) in terms of Random Matrix Theory (RMT) 
and Random Plane Wave (RPW) modeling. We are therefore 
in the position to follow the line of reasoning in Ref. 13 to 
analyse the SDFT calculation. We shall do this in this sec- 
tion first to get some understanding of the peak spacing and 
spin distributions obtained, and in a later stage to address the 
mechanism of spin contamination. 



A. Universal Hamiltonian form 

To model the statistical properties obtained from the SDFT 
calculations, let us assume that the Thomas-Fermi density 
across the dot has variation small enough that we can take the 
parameter 7-3 as a constant. We furthermore impose n'^-p (r) = 



nf^p(r) = nTF(r)/2 ( J n^p(r)dr might then be half integer 
for odd N, which is not a problem since at the Thomas-Fermi 
level the quantization of particle number is not playing any 
role) and write the second order Strutinsky correction as 



(45) 



with /io- = 0, 1 the occupation number of orbital i with spin 

(46) 



and 



Mff = I dvdv' \Uv)\^ 



(7.17 

scrccnci 



,(r,r')|0,(r')|^ (47) 



For a chaotic system, it can be shown that the fluctuations 
of the M"'^ are small (variance of order ^ \ng/ g^) and that 
their mean values are given by 



(T,CJ 

i,3 



scrccnc( 



,(g-o) + MT5y(v; 



(J,CJ \ 

screened/ 



lA 



(48) 

where A is the area of the dot, /it is 2 here since time reversal 
symmetry is preserved (but would be 1 if it were broken), and 



^ Screened (^^^ 



2ir 



= — / do v: 



1 

2n 



screened 



{^/2{l + cos 9)kF) 



(49) 

is the average of the screened interaction over the Fermi circle 
(note a^kp = \/2/rs). 

Note however that the screened interaction Eq. ( I26> is de- 
rived under the assumption that the oscillating part of the den- 
sity integrates to zero, so that the total displaced charge pro- 
viding the screening also does. Between the reference S — 
configuration (the TF case) and the higher S ones, the total 
number of electrons is, of course, conserved, but not the num- 
ber for each spin. It should be kept in mind, therefore, that the 
q = component of the density cannot be screened; this can 
be included simply by setting \s{q — 0) = X]^^'"^{q = 0). 

The screened interaction matrix is characterized by its 
eigenvalues Eqs. j40> and J41> . In Fig.|5] we thus plot the 
dependence of these quantities averaged over the Fermi circle, 
and compare them to their bare counterparts. 

A few remarks are in order concerning this figure. First, 
we see that, because of the divergence of the Coulomb inter- 
action at small q, screening has a drastic effect for the charge 
channel. Screening is less dramatic in the spin channel, but 
can still be an order one effect as increases. 

Furthermore, while the screening decreases the absolute 
strength of the interaction in the charge channel, it actually in- 
creases it in the spin channel. Indeed, since the interaction in 
the spin channel coming from SDFT is attractive, the charges 
in the bulk of the Fermi sea will, as long as this does not in- 
volve too much kinetic energy, move so as to increase the spin 
polarization. 
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FIG. 9: (Color online) Average over the Fermi circle of the eigenval- 
ues of the screened and bare SDFT interactions [in units of 2/7V(0)] 
as a function of = {%a^n)~^^'^ . Dark: charge channel. Lighter 
color (green online): spin channel. Dashed: bare interaction Q^'^s'^) 
(with a cutoff of the momentum at g ~ 1/L in the charge channel). 
Solid: screened interaction (A^s) with rj — 0.25. Thin dot-dashed: 
same but with r; = 0. Since the q dependence of As is entirely due to 
the Tr^p correction to the Thomas-Fermi kinetic energy functional, it 
therefore disappears in this latter case. The interaction in the charge 
channel is, of course, dramatically decreased by screening; in con- 
trast, screening increases the magnitude of the interaction in the spin 
chan 
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FIG. 10: (Color online) Mean value of the gain in interaction energy 
- (AE'-^'>{S, + 1) - A£'(^)(5'^)) (thick lines) and of the mean 
one-particle energy cost (thin horizontal lines) associated with flip- 
ping the spin of one particle in the quantum dot. Solid: Sz = 0; 
Dot-Dashed = 1/2; Dashed: = 1 



Finally, we note that for the value of the parameter rj that we 
use, the effect of the first h correction T^p on the screened in- 
teraction is very small in the charge channel, and only slightly 
larger in the spin channel. 

If we neglect the fluctuations of the M^'j^ , Eqs. ( I39> , ( I45> . 
and ( I48> imply that Ai?^^^ is just a function of the number 
of particles N — Na + Np in the dot and the z-component 
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FIG. 11: (Color online) Comparison between the analytical RPW 
predictions [Eqs. <48> and <51H and numerical calculations of the 
mean and variance of Mij 's. The wave functions used in the numer- 
ical calculations are eigenfunctions of the effective Thomas-Fermi 
potential with A'^ — 200 in the quartic oscillator system. Lines 
(points) correspond to analytical (numerical) results, with the solid 
(circle) for M^f, short-dashed (square) for M°'f, long-dashed (up- 
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a, 13 



triangle) for M!^f and dot-dashed (down-triangle) for M°'l 



Sz — {Na — Np)/2 of the groundstate the spin. 



~ ^((Ac)fc — (As)fc)<5'z 



(50) 



The main value of this expression is how it compares to the 
universal Hamiltonian form,^^ -^ and we shall come back to 
this point in the discussion section. Already we can see, 
however, that it contains almost all the information neces- 
sary to understand qualitatively the ground-state spin distribu- 
tions. Indeed, looking at Fig.^| which shows the difference 
- (A£;(2) {N,Sz + l)~ A£:(2) (TV, Sz)) as a function of r, 
for several values of Sz, we see that for ~ 0.85, the in- 
teraction energy gain and one particle energy cost of forming 
a triplet are equal on average, and therefore triplets should 
become as probable as singlets. In the same way, spin 3/2 be- 
comes as probable as 1/2 at rg — 1.8, and spin 2 as probable 
as 1 at rs ~ 2.8. 

To have a more precise idea of the whether the random- 
plane-wave model captures the main physics, we can follow 
the approach of Ref. 13 and make a simulation of the peak 
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FIG. 12: (Color online) Peak spacing distributions for the 
RMT/RPW model corresponding to = 0.3 (left) and 1.3 (right). 
Solid : even TV, dashed : odd A'^. Inset : corresponding spin distribu- 
tion. 



spacing and spin distributions. We use GOE distributed en- 
ergy levels for the first-order correction Eq. i20\ . and take 
the second-order correction in the form Eq. ( I45> with the 

Af j j independent variables with mean Eq. (|48} and a vari- 
ance which can be computed using the method of Appendix A 
ofRef.H 



yaT[M^f] 32 



'''F^ dx 1 



[x) (51) 



IS? TTlkpLY X 4 

(x) + 5,, {v""' {x) + v""' (0) + (v/4-a;2) 



where 7)'^'^ (g/fc^) = ^(0)Kcrec„cd('?)/2- That the random- 
plane-wave model correctly describe the wave function statis- 
tics is illustrated in Fig. where the analytic expressions 
Eqs. ( I48t and ( 15 1> are compared, as a function of r^, to the 
result from the actual eigenfunctions derived from Eq. ( I16> . 
The analytic expression for the mean is expected to be quite 
reliable and hence the good agreement. The variance ( 15 1> is 
less accurate - because of the cutoff used, for instance - and 
so we consider the agreement in Fig.^Jquite good. 

Fig-Eldisplays the peak spacing and spin distribution for 
Ts — 0.3 (corresponding to [al]) and rs = 1.3 ([a4]) coming 
from a simulation in which the fluctuations of level spacing 
and of the M"'^ are included. We see that the qualitative be- 
havior observed in Fig |8]- and in particular the lower panel - 
is very well reproduced. Thus the RMT/RPW approach using 
the LSDA interaction is successful in comparison with the full 
SDFT calculation. 



cannot involve any spin contamination. Spin contamination is 
a way, in the SDFT calculations, to lower the energy of the 
system without changing the total z component of the spin Sz 
by having different wavefunctions for the a and /3 orbitals. 
However, the eigenstates implicit in the Strutinsky approach 
are almost identical to those of TJtf, 0f , and the are nearly 
independent of the spin. 

In this section, we shall discuss how this spin contamina- 
tion mechanism could be understood within this Strutinsky 
scheme. Rather than trying to consider the most general sit- 
uation, we will limit ourselves to the case of even number of 
particle N , and a ground state z component spin equal to zero. 

Let iJcff [jt-tf] = -^V^ + t/|F(r) be the Thomas-Fermi 
Hamiltonian defining the orbitals [see Eqs. (I14t and ( I16H . 
What we have done is to construct a approximate solution of 
the SDFT equations Eq. (|5|l as (r) = 'Y^=\ 4>j (r) plus 
some screening charge. In this respect, an important point 
was that the resulting potential change SW^ given by Eq. ( 13 H 
was such that the matrix element ((/)| |(5C/'^|(/)j) for i ^ j was 
negligible, to first order in 1/g. 

However, finding an approximate solution of the Kohn- 
Sham equation implies only that one has an extremum of the 
spin density functional, but not necessarily an absolute mini- 
mum. As pointed out earlier, modifications of the wavefunc- 
tions change the electronic density only if they mix occupied 
and unoccupied orbitals. Therefore, when searching for a new 
extremum of the spin density functional, with some chance to 
actually get the true minimum, a natural choice is to mix the 
last occupied orbital with the first unoccupied one (fif.j^i 
(for M = N/2 = Na = Np). 

Let us therefore look for approximations (pf to the KS 
wavefunctions defined by (pf ~ (j>f for i < M and 



Vli+i = -sin' 



i-Af +sinr0M+i (52) 

'^(j)M + C0s9'^(j)M+l , 



with possibly a different value of the angle 6"^ for the two 
spins a = a, p. For these wavefunctions, the Thomas -Fermi 
Hamiltonian has a matrix element 

{pIj\Htf\^Ij^i) = cosr sinr (fTM+i - hi) (53) 

in terms of the Thomas-Fermi energies e^vz+i and e^/. 

This change in the wave functions produces a modification 
of the densities Sn'^ = If^s^P — IV'P which, once screening is 
taken into account, itself generates a perturbation potential 



B. Spin contamination 



E 

c7'—a.^(3 



dv'V° 



3,(r,r')<5n-'(r') 



(54) 



As we saw in the previous sections, the statistical proper- 
ties of the model quantum dots obtained from the full SDFT 
computation are, at least up to ~ 1.3, well reproduced 
by the various forms of the Strutinsky approximation. How- 
ever, we also saw that spin contamination, when present in the 
SDFT calculations, was actually degrading the accuracy of the 
Strutinsky approximation on a case by case basis. Indeed, by 
construction, the Strutinsky approximation as we presented it 



(The modification of the other wavefunctions, except for 
screening this term, does not play a role here.) The 
self-consistent condition for the angles is that 

(^^,|i7TF|'PX/+i) + {^%\m%v)\^i,^,) = 

In order to find matrix elements of 5U'^ , let us consider for 
a moment the case a = a. Eq. \52l implies that 



3(r - 6P)p,l, + sin(r - 0^)^X,+i (55) 
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and therefore 



= 2cos(r 



I sm 



)/ (56) 



with 



/ = /drdrV?^(r)^^,+i(r)l/:;f_,(r,r')^^,(r'Xf+i(r') 



\ s 



a,/3 

screened 



'id A 



(57) 



where the last equality applies in lowest order in 1 /g, that is, 
neglecting the fluctuations. 

Noting finally that the same reasoning can be applied to the 
/3 orbitals, we find that to get an extremum for the spin density 
functional the angles 9°" and 9^ should obey 



cos 6*" sin6'"(eAf+i - e^/) = 
2cos(6'" -6''3)sin(6'" 



(58) 



and the analogous equation with a and jS interchanged. Obvi- 
ous solutions ai-e {9'^ = 9P = 0), (61" = 0, 6"^ = 7r/2), (6*" = 7r/2, 
9^ = 0), and (9°' = 9^ = 7r/2). The first one corresponds 
to the standard 5 = (non-spin-contaminated) solution; the 
other ones involve promotion of particles from the last occu- 
pied orbital to the first unoccupied one, and so are obviously 
of higher energy. However, other solutions may exist. Clearly 
they should fulfill 



008 6*" sinr = cos6l'^sin6''5 ; 
that is, up to an irrelevant multiple of tt phase 



9^ 



and therefore 



cos2r = 



£m+i — Cm 



(59) 



(60) 



(61) 



It can be checked that whenever the condition Eq. ( 16 U can 



raj) 



the 



be fulfilled [ie. when cm+i ~ < 4(V;°^^^,j^j;tcJ 
corresponding extremum has an energy smaller by a quantity 
sin'^(26'")[(As)fc — {Xc)fc]/A with respect to the non contam- 
inated configuration. In this situation, the spin contaminated 
Sz = state will be favored (but its energy still needs to be 
compared to the lowest energy state with Sz = 1). 



VI. DISCUSSION 

To summarize our findings, we have seen that up to val- 
ues of order one (in practice, 1.3 here), the Strutinsky approxi- 
mation yields a ground state total energy with fluctuating part 
reliable up to typically five percent of the mean level spac- 
ing. Furthermore, these errors can be related to the occur- 
rence of spin contamination in the SDFT calculation, which 
cannot be reproduced by the Strutinsky scheme (in its sim- 
plest form). Indeed, as discussed in the last section, this latter 
gives, by construction, an approximation to an extremum of 
the Kohn-Sham functional for which the a and [3 orbitals are 




1 
Spacing 

FIG. 13: (Color online) Same as Fig. ll2r b). but with a residual inter- 
action modeled by a perturbative calculation in the interaction poten- 
tial VRPA(r — r') whose parameters are set by Eqs. <64t and <65> . 



nearly identical, but not necessarily an approximation to the 
true minimum. 

Nevertheless, the qualitative properties of the peak spacing 
and spin distribution are correctly reproduced by the Strutin- 
sky approximation (at least up to Vg of order one, for which 
spin contamination does not appear to change drastically the 
distributions). For a chaotic confining potential this makes it 
possible to use the modeling in terms of Random Matrix The- 
ory (for the energy levels) and Random Plane Waves (for the 
eigenstates) introduced in Ref. 13. Within this RMT/RPW 
modeling, and in the limit of large dots for which fluctuations 
of the residual interaction term are small, the main features of 
the peak spacing and spin distributions can be understood as 
arising from the interplay between one particle level fluctua- 
tions and the spin dependence of the mean residual interaction 
term [see Eq. ( 15 OH 



(2) 

SDFT 



{N,Sz 



\ bare q2 
Ac O ^ 



((Ac)fc - {Xs)tc)Sz 



(term depending on N only) . (62) 



It is useful to compare this expression, as well as the cor- 
responding distributions based on RMT/RPW modeling such 
as those in Fig. [21 (which actually take into account the fluc- 
tuations of the residual interaction term), to what is obtained 
following the more traditional route to the analysis of peak 
spacing and spin distributions for chaotic quantum dots.iSii^ii^ 
In these earlier approaches, the ground state energy of the 
quantum dot would, in a way very similar way to Eq. ( I19> . be 
described as the sum of a large non-fluctuating classical-like 
term, a one-particle energy contribution computed for some 
effective confining potential, and finally a residual interaction 
term AiJ^p^. This latter would, however, be understood as 
originating from a weak interaction VRpA(r — r') which in 
the Random Phase Approximation can be shown to be just the 
RPA screened potential|i2ii4 but should be understood more 
properly as the residual interaction between quasi-particles in 
Landau Fermi liquid spirit. We shall thus refer to this latter 
description as the "RPA" approach, although this is slightly 
inappropriate. 

While the Strutinsky approximation to SDFT gives a resid- 
ual interaction which can be understood as a first-order per- 
turbation (without exchange) in terms of the spin dependent 
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model 


S = S^i S = l 5'=| 5 = 2 


SOFT 
ST* 
ST*/RPW 
RPA/RPW (unscreened Cooper) 
RPA/RPW (screened Cooper ) 


0.42 ±0.03 0.76 ±0.03 0.54 ± 0.03 0.23 ± 0.03 0.03 ±0.01 
0.34 ± 0.03 0.74 ± 0.03 0.61 ± 0.03 0.25± 0.03 0.04 ± 0.01 
0.28 0.68 0.62 0.30 0.09 
0.48 0.84 0.49 0.16 0.03 
0.58 0.89 0.40 0.11 0.02 



TABLE I: Ground state spin probability for various model introduced in this paper. The two first lines correspond respectively to the SDFT 
calculation and [second form of] Strutinsky approximation for the quartic oscillator systems introduced in section lTvl The statistics are built 
from the ground state spin of dots containing 100 to 200 electrons, for a few confining potential corresponding to an interaction parameter 
Vs — 1.3. The three last lines are the results of RMT/RPW modeling for 150 electron and the same value of the interaction parameter rs, 
and respectively: the interaction derived from SDFT (third line); the "RPA"-like interaction using a screened Cooper channel, i.e. such that 
Eqs. <64t and <65> apply (fourth line); and the "RPA"-like interaction using a unscreened Cooper channel, i.e. such that Eqs. <64t applies but 
C = Js (fifth line). 



potential Eq. i26\ . in contrast one has in the "RPA" approach 
a residual interaction arising from the perturbative corrections 
in some VRPA(r, r') (including both direct and exchange) as 
well as possibly higher-order terms which turn out to be im- 
portant for time-reversal invariant systems (the Cooper series). 
Under this assumption and following exactly the same analy- 
sis as the one leading to Eq. ( I62l i. one gets for the mean resid- 
ual interaction term in the RPA approach 



AE^^l^iN, S,) = JsS,{S, + 1) ± C(mt - l)S. 

+ (term depending on N only) (63) 

(again /ix = 2 but would be 1 if time-reversal invariance were 
broken). 

In Eq. (I63> . the parameter Js is equal to —{V{q)){c, where 
the Fermi circle average is defined by Eq. i49i and VRPA(q) 
is the Fourier transform of Vrpa (r — r' ) . More properly how- 
ever, one should understand Js as being related to Fermi liq- 
uid parameter /q*^' through 

Js/A = fl,"^ . (64) 

In first-order perturbation theory, ( would be equal to Jg, but 
screening associated with higher oder terms in the Cooper 
channel somewhat reduces this value. An analysis following 
the lines of Ref. 30 suggests that 

C = rrr^ ' (65) 

with kp the Fermi momentum and L the typical size of the 
system. We shall assume this in the following discussion, 
bearing, in mind that this is true only "up to logarithmic accu- 
racy". 

The remarkable point here is that X^^'-'''{rs)/A is actually 

the same thing as /g"' (rs ), in the sense that the Landau Fermi 

liquid parameter ,fQ°'\rs) can be interpreted as the second 
derivative with respect to the polarization, at fixed total den- 
sity, of fixe the exchange correlation energy per particle of the 
uniform electron gas. This implies that the term quadratic in 
Sz of Eqs. i62\ and ( I63t actually do correspond. 



As a consequence if we compare Fig. I12f b). which is ob- 
tained from a RMT/RPW simulation with the interaction cor- 
responding to the spin density functional, with Fig. ob- 
tained in the same way but with an "RPA"-like interactionpi 
the differences in spin polarization and in odd/even asymme- 
try for the peak spacing distribution can be almost entirely 
associated with the different linear Sz terms in Eqs. (|62} and 

To get some further insight into the difference between the 
two approaches, let us consider Table|I]which shows the spin 
distributions for an interaction parameter r<; ~ 1.3 and various 
approximations discussed in this paper. Comparing different 
pairs of lines gives a sense of the importance of the various 
issues. For instance, comparing the second line with the third 
gives an idea of how accurate the RMT/RPW model is for 
the statistical properties of the real energy levels and eigen- 
functions of the quartic oscillator system, since both lines are 
based on the same spin dependent interaction Eq. ( I26t . The 
first line compared to the second, on the other hand, is a mea- 
sure of how well the Strutinsky scheme approximates the full 
SDFT calculation. Presumably, most of the difference be- 
tween these two lines can be associated with the presence of 
spin contamination in SDFT - it is in some sense a measure 
of the effectiveness of spin contamination in lowering the total 
spin of the system. Going further down the table, the differ- 
ence between the third and the two last lines is a measure of 
the impact of different linear terms in Sz in Eqs. i62\ and i63i . 
without screening the Cooper channel for the fourth line and 
with a screened Cooper channel [according to Eq. ( I65l ll for 
the last line of the table. 

From Table U it appears that within the accuracy of the 
RMT/RPW modeling, which seems to be around 5%, the 
SDFT result is compatible with an RPA-like approach if the 
Cooper channel is not screened. In other words, the fact that 
((As)fc — (Ac)fc)/2 is more negative that Js = A^*^™ - pro- 
ducing higher spins - is compensated by the effect of spin 
contamination - which favors lower spins. As seen in Ref.l3ll 
this compensation between the two effects seems to exist also 
for higher values of r-g. On the other hand, spin contamina- 
tion is not a sufficiently strong effect to compensate for the 
absence of screening of the Cooper channel. 

It remains to decide which of the two approaches is the 
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more correct. This, in the end, can only be addressed by 
comparison with exact calculations for quantum dots (e.g. 
quantum Monte Carlo). One argument that may be consid- 
ered is that in the presence of a time-reversal breaking term 
(i.e. — 1), general symmetry considerations impose that 
the mean value of the residual interaction term is a function 
of S{S + 1), but not independently of and S. Expres- 
sions (I62> clearly do not fulfill this constraint, while Eq. (I63> 
does. Since, however, spin contamination seems to compen- 
sate for this difference it might just be that for time-reversal 
non-symmetric systems, SDFT and RPA basically agree. On 



the other hand, the screening of the Cooper channel does not 
seem to be reproduced by the SDFT calculations, and this 
might be the cause of the higher spin found in this approach. 
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For the higher density case [al], it was actually necessary to place 
the gate quite close to the electron gas to see any effect on the 
density, even on the classical scale. The image charge associated 
with the gate is at a distance 2zo which is, however, still larger 
than the screening length within the electron gas. 
Actually (Ac^'^°)fc would be divergent if we did not use a cutoff 
on q at the inverse of the system size. 

For this simulation, the fluctuation of the residual interaction term 
also depends on the precise way the Cooper channel contribution 
is screened. Fig. ll2l corresDonds to an evaluation of the variance 
of these fluctuations for which we have assumed no screening, 
and is therefore rather an upper bound. The fluctuations being in 
any case rather small, their magnitude does not change the peak 
spacing distribution drastically (their main effect is to broaden the 
sharp peak for odd spacings), and a more detailed treatment of the 
effect of screening of the Cooper channel on the fluctuations of 
the residual interaction term should not change significantly the 
picture. 



